      Subroutine Regroup

      include "param.fi"
      include "common.fi"

      If (DIMM .EQ. DIMM2D) Call TrackBall2D            

      ! Method 1: regroup two particles
      If (DIMM .EQ. DIMM2D) Call Regroup2D

      ! Method 2: Remeshing when any particle is bad
      !If (DIMM .EQ. DIMM2D) Call Remeshing2D

      End Subroutine

C=================================================================

      Subroutine Remeshing2D
      include "param.fi"
      include "common.fi"

      Integer I, J, K, OverLimitNum, II, JJ, IB, IB2
      Integer HI, HJ, HK, HII, HJJ, HKK
      Integer DELI, DELJ, DELK, IP
      Real ShapeCoef
      Real DirI, DirJ, Dir_Distance, tmp, tmp1
      Real NewX(NTOT3T, 3), NewV(NTOT3T, 3)
      Real NewP(NTOT3T), NewRho(NTOT3T)
      Real Center(3)

      ! Find the over-limited particles
      OverLimitNum = 0
      Do I = 1, NTOT3T
        ShapeCoef = ShapeP(I, 3) / ShapeP(I, 1)
        If (ShapeCoef .GT. ShapeLimit) Then
          OverLimitNum = OverLimitNum + 1
        Endif
      Enddo
      If (OverLimitNum .EQ. 0) Return

      Write(*,*) "Over-limited Particles:", OverLimitNum, "Remeshing..."

          Do IB = 1, NTOT3T
            Call NP2DIndex13(IB, I,J,K)
            ! New position
            NewX(IB, 1) = XMIN + DX * (REAL(I)-0.5)
            NewX(IB, 2) = 0.0
            NewX(IB, 3) = XMIN + DX * (REAL(K)-0.5)
            !Write(*,*) IB, I, K

            Call Interpolation(NewX(IB,1), NewX(IB,2), NewX(IB,3),
     &                         U1, NewV(IB,1))
            Call Interpolation(NewX(IB,1), NewX(IB,2), NewX(IB,3),
     &                         U2, NewV(IB,2))
            Call Interpolation(NewX(IB,1), NewX(IB,2), NewX(IB,3),
     &                         U3, NewV(IB,3))
            Call Interpolation(NewX(IB,1), NewX(IB,2), NewX(IB,3),
     &                         P, NewP(IB))
            Call Interpolation(NewX(IB,1), NewX(IB,2), NewX(IB,3),
     &                         Rho, NewRho(IB))
          Enddo ! IB
          Do I = 1, NTOT3T
            ! 3) Overwrite the old particles
            X1(I) = NewX(I, 1)
            X2(I) = NewX(I, 2)
            X3(I) = NewX(I, 3)
            U1(I) = NewV(I, 1)
            U2(I) = NewV(I, 2)
            U3(I) = NewV(I, 3)
            P(I)  = NewP(I)
            Rho(I)= NewRho(I)
            ! 4) Don't forget to reset the shape coefficients
            ShapeP(I,1) =0.5*Sqrt(Weight(I))
            ShapeP(I,3) =0.5*Sqrt(Weight(I))
            ShapeVec2D(I,1)=1
            ShapeVec2D(I,2)=0
            ShapeVec2D(I,3)=0
            ShapeVec2D(I,4)=1
          Enddo

      End Subroutine


C=================================================================

      Subroutine Regroup2D
      include "param.fi"
      include "common.fi"

      Integer I, J, OverLimitNum, II, JJ, IB, IB2, III, tag
      Integer HI, HJ, HK, HII, HJJ, HKK
      Integer DELI, DELJ, DELK, IP
      Integer OverLimitIndex(NTOT3T)
      Real ShapeCoef
      Integer JIndex(30), JNum, JDir(30), Jmin
      Real DirI, DirJ, Dir_Distance, tmp, tmp1
      Real NewX(2, 3), NewV(2, 3), NewP(2), NewRho(2)
      Real Center(3)
      Real Coef_Axis

      !Coef_Axis = 4.0
      Coef_Axis = 3.88 ! Correction from test cases

      ! Find the over-limited particles
      OverLimitNum = 0
      Do I = 1, NTOT3T
        ShapeCoef = ShapeP(I, 3) / ShapeP(I, 1)
        If (ShapeCoef .GT. ShapeLimit) Then
          OverLimitNum = OverLimitNum + 1
          OverLimitIndex(OverLimitNum) = I
          !Write(*,*) I, ShapeCoef, ShapeP(I, 3), ShapeP(I, 1)
        Endif
      Enddo
      If (OverLimitNum .EQ. 0) Return

      Write(*,*) "Over-limited Particles:", OverLimitNum

      ! For eachone, find similar particles
      Do I = 1, OverLimitNum
        II = OverLimitIndex(I)
        If (II .EQ. 0) Cycle
        DirI = ATAN(ShapeVec2D(II, 2) / ShapeVec2D(II, 1))
        JNum = 0
        Do J = I + 1, OverLimitNum
          JJ = OverLimitIndex(J)
          If (JJ .EQ. 0) Cycle
          ! Direction should be similar
          DirJ = ATAN(ShapeVec2D(JJ, 2) / ShapeVec2D(JJ, 1))
          Dir_Distance = ABS(DirJ - DirI)
          If (Dir_Distance .GT. 3.14/2.0) Then
            Dir_Distance = 3.14 - Dir_Distance
          Endif
          If (Dir_Distance .GT. 0.2) Cycle
          ! They should be near to each other
          tmp = Distance2P(X1(II), X2(II), X3(II),
     &                     X1(JJ), X2(JJ), X3(JJ))
          If (tmp .GT. ShapeP(II, 1) * 3.6) Cycle
          ! The long axis should be similar
          If (ShapeP(II,3).GT.ShapeP(JJ,3)*1.3 .OR.
     &        ShapeP(JJ,3).GT.ShapeP(II,3)*1.3) Cycle
          ! There should be no particle near the new positions
          Center(1) = PeriodMid(X1(II), X1(JJ))
          Center(3) = PeriodMid(X3(II), X3(JJ))
          NewX(1, 1) = InBox(Center(1) + ShapeVec2D(II, 3) *
     &        (ShapeP(II,3) + ShapeP(JJ,3)) / Coef_Axis)
          NewX(1, 2) = 0.0
          NewX(1, 3) = InBox(Center(3) + ShapeVec2D(II, 4) *
     &        (ShapeP(II,3) + ShapeP(JJ,3)) / Coef_Axis)
          NewX(2, 1) = InBox(Center(1) - ShapeVec2D(II, 3) *
     &        (ShapeP(II,3) + ShapeP(JJ,3)) / Coef_Axis)
          NewX(2, 2) = 0.0
          NewX(2, 3) = InBox(Center(3) - ShapeVec2D(II, 4) *
     &        (ShapeP(II,3) + ShapeP(JJ,3)) / Coef_Axis)
          tmp = MIN( 0.25 * (ShapeP(II,3) + ShapeP(JJ,3)),
     &               1.0 * (ShapeP(II,1) + ShapeP(JJ,1)) )
          tag = 0
          Do III = 1, NTOT3T
            If (III.NE.II.AND.III.NE.JJ) Then
              tmp1 = Distance2P(NewX(1,1), NewX(1,2), NewX(1,3),
     &                          X1(III),   X2(III),   X3(III))
              If (tmp1.LT.tmp) Then
                tag = 1
                Exit
              Endif
              tmp1 = Distance2P(NewX(2,1), NewX(2,2), NewX(2,3),
     &                          X1(III),   X2(III),   X3(III))
              If (tmp1.LT.tmp) Then
                tag = 1
                Exit
              Endif
            Endif
          Enddo
          If (tag.EQ.1) Cycle
          ! Add it to the List
          JNum = JNum + 1
          JDir(JNum) = Dir_Distance
          JIndex(JNum) = J
        Enddo
        ! If there are several particles, find the best one
        If (JNum .GT. 0) Then
          ! Find the one with smallest Dir_Distance
          tmp = 10.0
          Do J = 1, JNum
            If (JDir(J) .LT. tmp) Then
              tmp = JDir(J)
              Jmin = J
            Endif
          Enddo
          ! There must be a value in Jmin
          II = OverLimitIndex(I)
          JJ = OverLimitIndex(JIndex(Jmin))
 
!          Write(*,*) "=========== Change", II, JJ
!          Write(*,*) X1(II), X3(II), Rho(II), P(II), U1(II),
!     &               U3(II), Weight(II)
!          Write(*,*) X1(JJ), X3(JJ), Rho(JJ), P(JJ), U1(JJ),
!     &               U3(JJ), Weight(JJ)

          ! To regroup them
          ! In OverLimitNum: I, JIndex(Jmin)
          ! Overall Index  : OverLimitIndex(I),
          !                  OverLimitIndex(JIndex(Jmin))
          !         =====>   II and JJ
          ! 1) Calculate the locations of new particles
          Center(1) = PeriodMid(X1(II), X1(JJ))
          Center(3) = PeriodMid(X3(II), X3(JJ))
          NewX(1, 1) = InBox(Center(1) + ShapeVec2D(II, 3) *
     &        (ShapeP(II,3) + ShapeP(JJ,3)) / Coef_Axis)
          NewX(1, 2) = 0.0
          NewX(1, 3) = InBox(Center(3) + ShapeVec2D(II, 4) *
     &        (ShapeP(II,3) + ShapeP(JJ,3)) / Coef_Axis)
          NewX(2, 1) = InBox(Center(1) - ShapeVec2D(II, 3) *
     &        (ShapeP(II,3) + ShapeP(JJ,3)) / Coef_Axis)
          NewX(2, 2) = 0.0
          NewX(2, 3) = InBox(Center(3) - ShapeVec2D(II, 4) *
     &        (ShapeP(II,3) + ShapeP(JJ,3)) / Coef_Axis)
          tmp = MIN(0.25 * (ShapeP(II,3) + ShapeP(JJ,3))
     &             , 1.0 * (ShapeP(II,1) + ShapeP(JJ,1)))
          ! 2) Make interpolation on the new particles
          Do IB = 1, 2
            !! Simplest method: average
            !NewV(IB, 1) = 0.5 * (U1(II) + U1(JJ))
            !NewV(IB, 2) = 0.5 * (U2(II) + U2(JJ))
            !NewV(IB, 3) = 0.5 * (U3(II) + U3(JJ))
            !NewP(IB)    = 0.5 * (P(II) + P(JJ))
            !NewRho(IB)  = 0.5 * (Rho(II) + Rho(JJ))
            Call Interpolation(NewX(IB,1), NewX(IB,2), NewX(IB,3),
     &                         U1, NewV(IB,1))
            Call Interpolation(NewX(IB,1), NewX(IB,2), NewX(IB,3),
     &                         U2, NewV(IB,2))
            Call Interpolation(NewX(IB,1), NewX(IB,2), NewX(IB,3),
     &                         U3, NewV(IB,3))
            Call Interpolation(NewX(IB,1), NewX(IB,2), NewX(IB,3),
     &                         P, NewP(IB))
            Call Interpolation(NewX(IB,1), NewX(IB,2), NewX(IB,3),
     &                         Rho, NewRho(IB))
          Enddo ! IB
!          Write(10,*) ShapeP(II,1), ShapeP(II,3), 
!     &                ShapeP(JJ,1), ShapeP(JJ,3)
          ! 3) Overwrite the old particles
          X1(II) = NewX(1, 1)
          X2(II) = NewX(1, 2)
          X3(II) = NewX(1, 3)
          U1(II) = NewV(1, 1)
          U2(II) = NewV(1, 2)
          U3(II) = NewV(1, 3)
          P(II)  = NewP(1)
          Rho(II)= NewRho(1)
          X1(JJ) = NewX(2, 1)
          X2(JJ) = NewX(2, 2)
          X3(JJ) = NewX(2, 3)
          U1(JJ) = NewV(2, 1)
          U2(JJ) = NewV(2, 2)
          U3(JJ) = NewV(2, 3)
          P(JJ)  = NewP(2)
          Rho(JJ)= NewRho(2)
          ! 4) Don't forget to reset the shape coefficients
          tmp = 0.25 * (ShapeP(II,3) + ShapeP(JJ,3))
          tmp1 = 1.0 * (ShapeP(II,1) + ShapeP(JJ,1))
          If (tmp1>tmp) Then 
	    ShapeP(II, 3) = tmp
          ShapeP(II, 1) = tmp1
          ShapeP(JJ, 3) = tmp
          ShapeP(JJ, 1) = tmp1
	    Else 
	    ShapeP(II, 3) = tmp1
          ShapeP(II, 1) = tmp
          ShapeP(JJ, 3) = tmp1
          ShapeP(JJ, 1) = tmp
	    Endif

          !Write(*,*) Rho(II), P(II), U1(II)
!          Write(*,*) "=========== After Change"
!          Write(*,*) X1(II), X3(II), Rho(II), P(II), U1(II),
!     &               U3(II), Weight(II)
!          Write(*,*) X1(JJ), X3(JJ), Rho(JJ), P(JJ), U1(JJ),
!     &               U3(JJ), Weight(JJ)


          ! Mark them
          OverLimitIndex(I) = 0
          OverLimitIndex(JIndex(Jmin)) = 0
        Endif
      Enddo

      ! see how many particles left
      J = 0
      Do I = 1, OverLimitNum
         If (OverLimitIndex(I) .GT. 0) Then
           J = J + 1
         Endif
      Enddo
      Write(*,*) "After regroup:", J

      !Pause
      

      End Subroutine

C=================================================================

      Subroutine TrackBall2D
      ! ShapeP(I, 1) -- x length of the particle
      ! ShapeP(I, 3) -- z length of the particle

      include "param.fi"
      include "common.fi"

      Integer I
      Real PRadius
     
      real :: S(2,2),A(2,2)

      type :: Oval
        real :: l(2,1),r(2,1)
        real :: a2,b2
      end type  Oval
      type(Oval) :: OldE,NewE

      If (ifDual.EQ.1) Then
        Do I = 1, GridDual2
          Call Interpolation(X1Dual(I), 0.0, X3Dual(I),
     &                       U1, U1Dual(I))
          Call Interpolation(X1Dual(I), 0.0, X3Dual(I),
     &                       U3, U3Dual(I))
        Enddo
        Do I = 1, NTOT3T
          Call IntepGrad(X1(I), X3(I), U1Dual, DU1DX(I), DU1DZ(I))
          Call IntepGrad(X1(I), X3(I), U3Dual, DU3DX(I), DU3DZ(I))
        Enddo
      Else
        Call Grad(U1, DU1DX, DU1DY, DU1DZ)   ! DU1DY is dummy
        Call Grad(U3, DU3DX, DU3DY, DU3DZ)   ! DU3DY is dummy
      Endif

      if(tagRegroup.EQ.1)then
      !----------------------------
        Do I = 1, NTOT3T
          ShapeP(I, 1) = ShapeP(I, 1) + DU1DX(I) * DT * ShapeP(I, 1)
          ShapeP(I, 3) = ShapeP(I, 3) + DU3DZ(I) * DT * ShapeP(I, 3)
        Enddo
      !-----------------------------------
      else if(tagRegroup.EQ.2)then
        Do I=1,NTOT3T
          !--S----------------------
          S(1,1)=DU1DX(I)
          S(1,2)=0.5*(DU1DZ(I)+DU3DX(I))
          S(2,1)=0.5*(DU1DZ(I)+DU3DX(I))
          S(2,2)=DU3DZ(I)
          !--A---------------------
          A(1,1)=0
          A(1,2)=0.5*(DU1DZ(I)-DU3DX(I))
          A(2,1)=-0.5*(DU1DZ(I)-DU3DX(I))
          A(2,2)=0
          !-----------------------
          OldE%a2=1/ShapeP(I,1)**2
          OldE%b2=1/ShapeP(I,3)**2
          OldE%l(1,1)=ShapeVec2D(I,1)
          OldE%l(2,1)=ShapeVec2D(I,2)
          OldE%r(1,1)=ShapeVec2D(I,3)
          OldE%r(2,1)=ShapeVec2D(I,4)
          !----------------------------
          call ellipse(S,A,OldE,NewE,DT)
          ShapeP(I,1)=Sqrt(1/NewE%a2)
          ShapeP(I,3)=Sqrt(1/NewE%b2)
          ShapeVec2D(I,1)=NewE%l(1,1)
          ShapeVec2D(I,2)=NewE%l(2,1)
          ShapeVec2D(I,3)=NewE%r(1,1)
          Shapevec2D(I,4)=NewE%r(2,1)
        Enddo
      endif
      End Subroutine

